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Abstract: A new filter named the maximum likelihood-based iterated divided difference 
filter (MLIDDF) is developed to improve the low state estimation accuracy of nonlinear 
state estimation due to large initial estimation errors and nonlinearity of measurement 
equations. The MLIDDF algorithm is derivative-fi"ee and implemented only by calculating 
the functional evaluations. The MLIDDF algorithm involves the use of the iteration 
measurement update and the current measurement, and the iteration termination criterion 
based on maximum likelihood is introduced in the measurement update step, so the 
MLIDDF is guaranteed to produce a sequence estimate that moves up the maximum 
likelihood surface. In a simulation, its performance is compared against that of the 
unscented Kalman filter (UKF), divided difference filter (DDF), iterated unscented Kalman 
filter (lUKF) and iterated divided difference filter (IDDF) both using a traditional iteration 
strategy. Simulation results demonstrate that the accumulated mean-square root error for 
the MLIDDF algorithm in position is reduced by 63% compared to that of UKF and DDF 
algorithms, and by 7% compared to that of lUKF and IDDF algorithms. The new algorithm 
thus has better state estimation accuracy and a fast convergence rate. 

Keywords: nonlinear state estimation; divided difference filter; maximum likelihood 
surface; target tracking 
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1. Introduction 

The problem of estimating the state of a nonhnear stochastic system from noisy measurement data has 
been the subject of considerable research interest during the past few years. Up to now the extended 
Kalman filter (EKF) has unquestionably been the dominating state estimation technique [1,2]. The EKF 
linearizes both the nonlinear process and the measurement dynamics with a first-order Taylor series 
expansion about the current state estimate. However, its accuracy depends heavily on the severity of 
nonlinearities. The EKF may introduce large errors and even give a divergent estimate when the 
nonlinearities become severe [3,4]. To improve the estimation accuracy, the second-order EKF 
proposed retains the Taylor series expansion up to the second term. The second-order EKF generally 
improves estimation accuracy, but at the expense of an increased computational burden [5]. Another 
attempt to improve the performance of the EKF involves the use of an iterative measurement update; 
the resulting algorithm is called the Iterated Extended Kalman filter (lEKF) [6]. The basic idea of 
lEKF is to linearize the measurement model around the updated state rather than the predicted state. 
This is achieved iteratively, and it involves the use of the current measurement. The lEKF has been 
proven to be more accurate on the condition that the state estimate is close enough to the true value, 
however, this is rarely the case in practice [7]. It was pointed out in [8] that the sequence of iterations 
generated by the lEKF and that generated by the Gauss-Newton method were identical, thus globally 
convergence was guaranteed. However, the Gauss-Newton method does not ensure that it goes up the 
likelihood surface [9,10]. Furthermore, EKF and lEKF require Jacobians, and the second-order KF 
requires Jacobians and Hessians. Calculation of Jacobians and Hessians is often numerically unstable 
and computationally intensive. In some system, the Jacobians and Hessians do not exit, which limits 
the applications of EKF, second-order EKF and lEKF. 

Recently, there has been development in derivative-free state estimators. The finite difference 
has been used in the Kalman filter framework and the resulting filter is referred to as the finite 
difference filter (FDF) [11]. The FDF uses the first-order difference to approximate the derivative of 
the nonlinear function; it may introduce large state estimation errors due to a high nonlinearity, similar 
to the EKF. The unscented Kalman filter (UKF) proposed in [12,13] uses a minimal set of 
deterministically chosen sample points to capture the mean and covariance of a Gaussian density. 
When propagated through a nonlinear function, these points capture the true mean and covariance up 
to a second-order of the nonlinear function. However, the parameters used in the UKF are required to 
tune finely in order to prevent the propagation of non-positive definite covariance matrix for a state 
vector's dimension higher than three. Another Gaussian filter, named the divided difference filter 
(DDF) was introduced in [14] using multidimensional Stirling's interpolation formula. It is shown 
in [15] that the UKF and DDF algorithms are commonly referred to as sigma point filters due to the 
properties of deterministic sampling and weighted statistical estimation [16], but the covariance 
obtained in the DDF is more accurate than that in the UKF. The iterated UKF with the variable step 
(lUKF-VS) in [10] proposed improved the accuracy of state estimation but its runtime was large due to 
its computation of the sigma points. Lastly, a relatively new technique called the particle filter (PF) 
uses a set of randomly chosen samples with associated weights to approximate the posterior 
density [17] and its variants are presented in [18]. The large number of samples required often makes 
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the use of PF computationally expensive, and the performance of PF is crucially dependent on the 
selection of the proposal distribution. Table 1 lists the pro and cons of the above filters. 

Table 1. Pro and cons for various filters. 



Various filters 



Advantage 



Disadvantage 



EKF 

Second-order EKF 
lEKF 
UKF 
lUKF-VS 
FDF 
DDF 
PF 



Less runtime 
high accuracy 
high accuracy 
Derivative-free 
Derivative-free, high accuracy 
Derivative-free, Less runtime 
Derivative-free 
Derivative-free 



Calculation of Jacobian 
Calculation of Jacobian and Hessian 
Calculation of Jacobian 
more runtime 
more runtime 
Lower accuracy 
Low accuracy 
Heavy computational burden 



The DDF also shows its weakness in the state estimation due to the large initial error and high 
nonlinearity in the application for state estimation of maneuvering target in the air-traffic control and 
ballistic re-entry target. Emboldened by the superiority of DDF, the basic idea of the lEKF and the 
iteration termination condition based on maximum likelihood, we propose a new filter named the 
maximum likelihood based iterated divided difference Kalman filter (MLIDDF). The performance of 
the state estimation for MLIDDF is greatly improved when involving the use of the iteration 
measurement update in the MLIDDF and the use of the current measurement. The remainder of this 
paper is organized as follows: in Section 2, we develop the maximum likelihood surface based iterated 
divided difference Kalman filter (MLIDDF). Section 3 presents the applications of the MLIDDF to 
state estimation for maneuvering targets in air-traffic control and ballistic target re-entry applications 
and discuss the simulation results. Finally, Section 4 concludes the paper and presents our outlook on 
future work. 

2. Development of Likelihood Surface Based Iterated Divided Difference Filter 

2.1. Divided Difference Filter 
Consider the nonlinear function: 

y = f(x) (1) 

Assuming that the random variable x £ R"^ has Gaussian density with mean x and covariance P^. 
The following linear transformation of x is introduced: 

z = S> (2) 

The transformation matrix S;c is selected as a square Cholesky factor of the covariance matrix such 
that = SxSx, so the elements of z become mutually uncorrelated [14]. And the function f is defined by: 

y = f(S,z) = f(z) (3) 

The multidimensional Stirling interpolation formula of Equation (3) about z up to second-order 
terms is given by: 
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y^~f{z) + bj + ^bli (4) 
The divided difference operators Di^^^, D^^f are defined as: 

^J = 7(ZAz^>"/.)f(^) (5) 

' j=i j=i ;=i ^ ^ 

i*] 

where Az = z - z and Azy is the y-th element of Az. / denotes a selected interval length, the optimal 
setting of / is selected as VS under the assumption that the estimation error are Gaussian and unbiased. 
The partial operators ^ and S are defined as: 



f(z + ^e^.) + f(z-^e^.) 



(7) 



^,f(z) = f(z + |e,)-f(z-|e,) (8) 



where ey is the unit column vector. 



We can obtain the approximate mean, covariance and cross-covariance of y using Equation (4): 

y = ^[y]-^f(x)+-^i:[f(x+/s,,.)+f(x-/s,,.)] (9) 
P^=^[(y-y)(y-yf] 

^' M (10) 



P^=£[(x-x)(y-yf]^i-Xs,,,(f(x + /s,,.)-f(x-/s,,,)f (n) 

^/ j=i 

where s^j is /-th column of Sx. 

Consider the state estimation problem of a nonlinear djoiamics system with additive noise, 
the %-dimensional state vector x,j of the system evolves according to the nonlinear stochastic 
difference equation: 

Xi=f(x,_i) + w,_i (12) 

and the measurement equation is given as: 

Zi=h(x,) + v, (13) 

Wfc_i and V/j are assumed i.i.d. and independent of current and past states, W/j;_i~J\r(0, Q/j;_i) , 
Vfe-JV^CCRfe). 
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Suppose the state distribution at k-\ time instant is X;j_i~J\r(Xfc_i, P^-i), and a square Cholesky 
factor of Pfc_i is Sx,k-i- The divided difference filter (DDF) obtained with Equations (9)-(l 1) can be 
described as follows: 

Step 1 . Time update 

(1) Calculate matrices containing the first- and second- divided difference on the estimated state 
Xfc_i at k-\ time: 

= {^[f (^.-1 + ^s,,,)-f(^.-. -^^,,)]} (14) 

Sg. = {^^[f(i,_, + /s,,) + f(i,_, -/s,,)-2f(x,_J]} (15) 

(2) Evaluate the predicted state and square root of corresponding co variance: 

= — ^f(Vi) + :^E[f(Vi+'s.j)+f(Vi-'Sxj)] (16) 

S,,,=rna([SS, S^,,_, Sgj) (17) 

Sy. j is 7-th column of S^^k-i ■ Tria{) is denoted as a general triagularization algorithm and Sy^ j^_-^ 
denotes a square-root factor of Q^-i such that Q^-i = ^w,k-i^w,k-i- 
Step 2. Measurement update 

(1) Calculate matrices containing the first- and second-divided difference on the predicted state x^: 

SS. = {^[h(x, +/\,)-h(x, -/s,,,)]} (18) 

SS, = {^^[h(x, +/s,,,) + h(x, (19) 

where Sxj is the j-th column of S;^ ^. 

(2) Evaluate the predicted measurement, square root of innovation covariance and cross-co variance: 

\ =^h(x,) + ^i:[h(x, +/s,,,)+h(x, -/s,,,)] (20) 
S^,,=rna([Sg, S,,, Sgj) (21) 

^.z,k='^.,ki^%ky (22) 

here Sv,/c denotes a square root of R/c such that = S^/jS^/j. 

(3) Evaluate the gain, state estimation and square root of corresponding covariance at k time: 

K,=(P^„,/SL,,)/S,,,, (23) 

^k-^k'^'^ki^k~'^k) (24) 
S.,.=rrm([[(S,^,-K,S«,) K,S,^, K,SL^!.]]) (25) 

here, the symbol "/" represents the matrix right divide operator. 
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2.2. Refining the Measurement Update Based on Divided Difference 



Consider X;^ and current measurement Zj^ as realization of independent random vectors with 
multivariate normal distributions, e.g., Xfc~J\r(x,j, P,j) and z,j~J\r(h(x,j), R,j). For convenience, the 
two vectors are formed to a single augmented one Z = [z^^ x^^ ]^ . According to the independent 
assumption, we have: 

Z~Ar(g(x,),Q) 



Here: 



g(x,) = [h(x),x,f ,Q: 



R, 0 
0 R 



(26) 



(27) 



The update measurement problem becomes the one that computing the optimal state estimation and 
corresponding covariance given Z, g and Q. 
Defining the objective function: 

f(^) = i||r(^)f (28) 

where r(^) = S(Z — g(^)) ,and r(^) is the second-order differentiable. Then the above update 
measurement becomes clearly a non-linear least squares problem. 

Assuming the /-iterate is x^'-*, we can obtain the following equation [8]: 



= +K« (z, -h(x«) + H« (x, -x«)) 



(29) 



where K 



(0 



We know the sequence of iterates generated by the lEKF and that generated by the Gauss-Newton 
method were identical, thus globally convergent was guaranteed. The initial state X/j is included in the 
measurement update, the value of Xj^ has a direct and large effect on the final state estimation. When 
the measurement model fully observes the state, the estimated state x^*^ is more approximate to the true 
state than the predicted state x,^ [7]. Substituting x^by x^'-* into Equation (29), hence, the following 
iterative formula is obtained: 



x« + K^"(z,-h(x^'0) 



(30) 



Compared to Equation (29), the Equation (30) is simpler, and the two equations are identical when 
there is a single iteration. 
Now we consider the gain: 



Kf = P,Hr(HrP.H«+R,)" 



(31) 



where the terms 



Hf^PfcH^') + R;, and 



P/jH^'"*^ are approximate innovation covariance and 



cross-covariance obtained by linearizing the measurement equation: 



P^1=H«P,(H«)^+R, 



)('■) _' 



(32) 
(33) 
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The terms of Equations (32) and (33) are achieved by expanding the measurement Equation (13) up 
to a first-order Taylor term so that the linearized error is introduced due to the high-order truncated 
terms. As for the highly nonlinear measurement equation, the accuracy for state estimation is 
decreased if the linearized error is only propagated in the Equation (30). To decrease the propagated 
error, we can recalculate Equations (18)-(22) to obtain the terms P^^ , P^^ the following way: 

^ = {|r[h(x« + /sS)-h(x« -/slj)]} (34) 

= {^[h(x« ) + h(x« _/§(0 )_2h(x«)]} (35) 

S2.=2na([slS« S,, S^^Jf]) (36) 
p(0 cCOr (37) 

-^1!*: ~ ^i^^LV (38) 

where S^''' is a square Cholesky factor of the covariance P^'^ 
Hence, we can obtain the following iterative formula: 

xr=x«+P^^(P^')"[z,-h(x«)] (39) 



2.3. Maximum Likelihood Based Iteration Termination Criterion 



^k 



< f is used 



In the measurement update step of the lEKF algorithm, the inequality 

as the criterion to terminate the iteration procedure, where £ is the predetermined threshold. The 
threshold e is crucial to successfully using the lEKF algorithm, but selecting a proper value of e is 
difficult [10]. The sequence of iterations generated according to the above termination condition has 
the property of global convergence; however, it is not guaranteed to move up the likelihood surface, so 
an iteration termination criterion based on maximum likelihood surface is introduced. 

Consider x^'-* and Z;^ as the realization of independent random vectors with multivariate normal 
distributions, i.e., x^^~M(xp,, P^'"*), Zfc~J\r(h(Xfe), R^^). The likelihood function of the two vectors Xj^ 
and Zfc is defined as: 

A(x,,yJ = con5^exp{-0.5[(x,-x«)^P^')(x,-x«) + (z,-h(xJ)"R,'(z^ (40) 

where P]^"'^ = (P]^*^)-^ 

Meanwhile, the likelihood surface is defined as follows: 

J(x, ) = (xl' > - X, )^ P^-" (if - X, ) + (z, - h(x, ))^ R,' (z, - h(x, )) (41) 

We know that the solution that maximizes the likelihood function is equivalent to minimizing the 
cost function /(x^). The optimal value ofJ(xj^) is difficult to obtain, but the following inequality holds: 

J(xr>)<J(xf) (42) 
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We say that /(x^'^''^'') is close to the maximum likelihood surface than /(x^'''), equivalently, x^'^"'^'' 
has a more accurate approximation than x^'-* to the minimum value of /(x/^) [10]. Extending 
Equation (42) and using P^''' = S^'''s['''^, we immediately obtain the following inequality: 



(xr>)^(s«s«o-'xr' +(zr^)^R.'zr^ < ii^y^i^ (43) 



where '^and and are defined as: 





The sequence generated is guaranteed to go up the likelihood surface using Equation (43) as the 
criterion to iteration termination. 

2.4. Maximum Likelihood Based Iterated Divided Difference Filter 

We have now arrived at the central issue of this paper, namely, the maximum likelihood based 
iterated divided difference filter. Enlightened by the development of lEKF and the superiority of DDF, 
we can derive the maximum likelihood based iterated divided difference filter (MLIDDF) which 
involves the use of the iteration measurement update and the current measurement. But in view of the 
potential problems exhibited by the lEKF, we shall refine the covariance and cross-covariance based 
on divided difference and use the termination criterion which guarantees the sequence obtained moves 
up the maximum likelihood surface. The MLIDDF is described as follows: 

Step 1 . Time update 

Evaluate the predicted state x^ and square Cholesky factor S/j. of the corresponding covariance 
using the Equations (14)-(17). 

Step 2. Measurement update 

Let x^^-* = X/j and S^""* = S/j. Suppose that the «-th iterates are x^'-* and S^'''. 

(1) Evaluate the first- and second-order difference matrices using Equation (34) and (35). 

(2) Evaluate the square root of innovation covariance and cross-covariance: 




(46) 



(47) 



(3) Evaluate the gain: 



(48) 




(49) 



(50) 



where s^'^ is the y-th column of S^'"*. 
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Step 3. If the following inequality holds: 



(51) 



The iteration returns to Step 2; otherwise continue to Step 4. x^'^^'', z^'^^-'are defined in the 
Equations (44) and (45). 

Step 4. If the inequality is not satisfied or if i is too large (/ > iVmax), and the ultimate state estimation 
and square root of corresponding covariance at k time instant are: 



(52) 

(53) 



The MLIDDF algorithm has the virtues of free-derivative and better numerical stability. The 
measurement update of MLIDDF algorithm is transformed to a nonlinear least-square problem; the 
optimum state estimation and covariance are solved using Gauss-Newton method, so the MLIDDF 
algorithm has the same global convergence as the Gauss-Newton method. Moreover, the iteration 
termination condition that makes the sequence move up the maximum likelihood surface is used in the 
measurement update process. 

3. Simulation and Analysis 

In this section, we reported the experimental results obtained by applying the MLIDDF to the 
nonlinear state estimation of a maneuvering target in an air-traffic control scenario and a ballistic target 
re-entry scenario. To demonstrate the performance of the MLIDDF algorithm we compared its 
performance against the UKF, DDF, and the iterated UKF (lUKF) and iterated DDF (IDDF), both 
using a traditional iteration strategy. 



3.1. Maneuvering Target Tracking in the Air-Traffic Control Scenario 

We consider a typical air-traffic control scenario, where an aircraft executes a maneuvering turn in a 
horizontal plane at a constant, but unknown turn rate Q.The kinematics of the turning motion can be 
modeled by the following nonlinear process equation [2,19]: 



1 


sinQr 


0 


1-cosQr 


0 


Q 


Q. 


0 


cosQr 


0 


-smQ.T 


0 


0 


1-cosQ.T 


1 


smQ.T 


0 


Q 


Q. 


0 


sin Q.T 


0 


cosQ.T 


0 


0 


0 


0 


0 


1 



(54) 



where the state of the aircraft x = [x x yy £1]; x and y denote positions, and x and y denote velocities 
in the x and j directions, respectively; Tis the time-interval between two consecutive measurements; 
The process nose Wfe_i~J\r(0, Q) with a nonsingular covariance: 
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Q = 



0 
0 



0 
0 



0 ^ 
0 



n = 



>V3 7^72' 
T'/2 



(55) 



J 



The parameters 91 and ^2 are related to process noise intensities. A radar is fixed at the origin of the 
place and equipped to measure the range, r and bearing, 0. Hence, the measurement equation is written: 

r T \ 



I 2 , 2 
,tan"'(j,/x,)_ 



+ v,. 



(56) 



where the measurement noise V/j~J\r(0, R) and R = diag([crr ffg]). 

The parameters used in this simulation were the same as those in [19]. To tracking the maneuvering 
aircraft we use the proposed MLIDDF algorithm and compare its performance against the DDF. We 
use the root-mean square error (RMSE) of the position, velocity and turn rate to compare the 
performances of two nonlinear filters. For a fair comparison, we make 250 independent Monte Carlo 
runs. The RMSE in position at time k is defined as: 



n ^ 



(57) 



where (Xj^''*, y^''') and (Xj^''', y^''') ^re the true and estimated positions at the z'-th Monte Carlo run. 
Similarly to RMSE in position, we may also write formulas of RMSE in velocity and turn rate. 

Owing to the fact that the filters is sensitive to initial state estimation, Figures 1-3 show the RMSEs 
in position, velocity and turn rate, respectively for DDF and MLIDDF in an interval of 50-100 s. As can 
be seen in Figures 1-3, the MLIDDF significantly outperforms the DDF. 

Figure 1. RMSE in position for DDF and MLIDDF. 




LU 



60 



70 80 
Time/s 



90 



100 



Sensors 2012, 12 



8922 



Figure 2. RMSE in velocity for DDF and MLIDDF. 
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Figure 3. RMSE in turn rate for DDF and MLIDDF. 
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In order to analyze the impact of iteration numbers on performance of the MLIDDF algorithm, the 
MLIDDFs with various iteration numbers are applied to position estimation of maneuvering target. 
Figure 4 shows the RMSEs of DDF and MLIDDFs with various numbers in position. From Figure 4, 
the RMSE in position for the MLIDDF with iteration number 2 begins to decrease largely comparable 
to that of DDF. The RMSE in position for MLIDDF with iteration number 5 significantly reduce, and 
the MLIDDF algorithm has very fast convergence rate. The RMSE of the position for MLIDDF 
decreases slowly when the iteration number is greater than 8; the reason is that the sequence generated 
has basically reached the maximum likelihood surface when the iteration termination condition is met. 
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3.2. State Estimation of Reentry Ballistic Target 
3.2.1. Simulation Scene 



The relative location of the ballistic target re-entry (BTR) and the radar are shown in Figure 5. The 
inertial coordinate system (ECI-CS) OxiyjZj shown in Figure 5 is a right-handed system with the 
origin O at the Earth's center, axis Oxj pointing in the vernal equinox direction, axis Ozj pointing 
in the direction of the North Pole A^. Its fundamental plane Ox^y; coincides with the Earth's 
equatorial plane. 

Figure 5. Geometry of radar and BTR. 



f 

i 


A' 










/ y 








X 




K — ^ 














TT"*" 



■V/ 



Assume that the radar is situated at the surface of the Earth and considering the orthogonal 
coordinate reference system named East-North-Up coordinates system (ENU-CS) O^xyz has its origin 
at the location of the radar. In this system, z is directed along the local vertical and x andy lie in a local 
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horizontal plane, with x pointing east and y pointing north. Assuming that the Earth is spherical and 
non-rotating and the forces acting on the target are only gravity and drag [20], we can derive the 
following state equation according to the kinematics of the ballistic re-entry object in the reentry phase 
in ENU-CS: 



=«>x^-i+Gv|/(x^-i) + w 



k-\ 



(58) 



where X;^ = [a:^ Xj^ y/j j/j Zj^ Pj^Y is the state of ballistic target re-entry, and: 



where 



0 = 



> 0 0 0 
0^00 

0 0^0 

0 0 0 1 



T 

1 



TOO 

0 r 0 
0 0 T 
0 0 0 



rV2" 

T 



^-1 ^k-\ 



k-\ 



'k-\ 



■^Hk-X 'k-\ 

PiK-x)y ■ _M(Zk-l+K) 
^k-l 



"^k-l 



'k-l 



V,, 



k-l 



hk-i = h-i 



R 



(59) 



(60) 



(61) 



(62) 



T is the time interval between radar measurements, y5 is the ballistic coefficient (kg/m^), // and Re 
are the Earth's gravitational constant and average Earth radius, respectively. Below 90 km at height, 
the air density p(h) is approximately modeled as an exponentially decaying function of height, 
pQi) = Cie~'^2fi (j^j^ C2 are constant, specifically, c\ = 1.227, C2 = 1.093 x 10""* for h < 9,144 m, and 
ci = 1.754, C2 = 1.49 X 10""* for h > 9,144 m) [21]. 

Process noise is assumed to be white noise with zero mean; its covariance is approximately 
modeled as [2]: 



0 
0 
0 



0 

0 
0 



0 
0 

0 



0 
0 
0 

92^ 



(63) 



where qi (in m^/s^) and q2 (in kg^/m'^s) are the intensity of noise. 
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According to relative geometry, the measurement equation in the ENU-CS is described as: 

Zi=h(x,) + v^ (64) 

here = {Rj^ Ej^ A^Y , and: 

Rk = ^l4+yl+4+VR (65) 

= tan-' zj ^xl+yl + (66) 

4 = tan-'>;,/x,+v^ (67) 
The measurement noise v,j is assumed to be white noise with zero mean and covariance: 

'H,=diag^al al crj]) (68) 

where a^^, a^, are the error standard deviations of range, elevation and azimuth, respectively. It is 
independent of the process noise W/j and initial state Xg. 

3.2.2. Numerical Results and Analysis 

The parameters used in simulation were: r= 0.1 s, = 5 m^/s^, q2 = 5 kg^/m'^s.The initial position 
and magnitude of the velocity: xo = 232 km, yo = 232 km, zo = 90 km and vo = 3,000 m/s, and the initial 
elevation and azimuth angle: Eq = 7n/6 and Ao = Tc/4.The ballistic coefficient was selected as 
J3 = 4,000 kg/m\ The error standard deviations of the measurements were selected as cr^ = 100 m, 
(Tg- = 0.017 rad, and cr^ = 0.017 rad. A threshold £"=10 was set in the lUKF and IDDF, and a maximum 
iteration number A^^ax = 8 was predetermined. 

From the above parameters given, we can obtain the initial true state: 

Xo=[232km -mim/s 232km -mim/s 90km -1500m/ s 4000kg /m^] 

and we select the corresponding covariance as: Pq = diag[l00^ 50^ 100^ 50^ 100^ 50^ 200^] . 

To compare the performance of the various filter algorithms, we also use RMSE in the position, 
velocity and ballistic coefficient. The position RMSE at k time of the baUistic target re-entry is 
defined as: 

RMSE^ik) = -4')' -4")') (69) 

where (^fc^y^'^z^ '') and (^^ ^y^'^^j^''') are the true and estimated position at the i-th Monte Carlo 
run, N is the Monte Carlo runs. Similarity to the RMSE in position we may also write formulas of the 
RMSE in velocity and ballistic coefficient. 

Figures 6-8 shows the position, velocity and ballistic coefficient RMSEs, respectively, for the 
various filters in an interval of 15-58 s. The initial state estimate Xq is chosen randomly from 
Xo~J\r(xQ, Pq) in each run. All the filters are initialized with the same condition in each run. We make 
100 independent Monte Carlo runs. 
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Figure 6. Position RMSEs for various filters. 
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Figure 7. Velocity RMSEs for various filters. 
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Figure 8. Ballistic coefficient RMSEs for various filters. 
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From Figure 6 we can see that the position RMSE of the MLIDDF is much less than the UKF and 
DDF because of the use of the current measurement in the step of the iteration measurement update of 
the MLIDDF, and is less than the lUKF and IDDF owing to involving the proposed iteration strategy, 
so the estimates provided by the MLIDDF are markedly better than those of the UKF and DDF 
algorithms, and are better than those of lUKF and IDDF algorithms. The MLIDDF also shows a 
significant improvement over the other filters in the estimation of the velocity, as evidenced by Figure 7. 

As to the estimation of the ballistic coefficient, in the Figure 8, we can see that there is no 
improvement in the RMSE in the initial interval of the observation period (t < 35 s) because there is no 
effective information about it, while in the remaining period (35 s < t < 58 s) the baUistic coefficient 
RMSE is decreased because the effective information about the balUstic coefficient fi-om the latest 
measurement is fully used. Especially, Figure 8 illustrates that toward the end of the trajectory the 
estimates provided by the MLIDDF are markedly better than those of the UKF, DDF, lUKF and 
IDDF algorithms. 

Meanwhile, we observe from Figures 6-8 that the UKF and DDF have almost the same 
performance in the problem and the performance of lUKF and IDDF algorithms are almost identical. 

For further comparison of performance for the various filters, the average accumulated position 
mean-square root error (AMSRE) is defined as follows: 



1 ^ Ir 1 M 



(70) 



where M is the total number of measurement data points. The formulas of the velocity and ballistic 
coefficient AMSREs can be written similar to the AMSREp. The position, velocity and ballistic 
coefficient AMSREs for the various filters are listed in Table 2. Table 3 lists the runtimes of the 
various filters. 

Table 2. AMSREs for various filters. 



Various algorithms 


AMSREp(m) 


AMSRE/m/s) 


AMSRE 0(kg/m^) 


UKF 


2521.684 


329.911 


155.735 


DDF 


2521.573 


329.903 


155.276 


lUKF 


1035.340 


259.173 


149.603 


IDDF 


1035.273 


260.771 


149.756 


MLIDDF 


968.746 


255.916 


144.953 



Table 3. Runtimes of various filters. 



Various algorithms 


Ruiitime(s) 


UKF 


1.0840 


DDF 


0.2888 


lUKF 


1.9918 


IDDF 


0.5133 


MLIDDF 


1.3074 



From Table 2, it is seen that the position AMSRE for the MLIDDF algorithm is reduced by 62% 
compared to the UKF and DDF, by 7% compared to lUKF and IDDF. The velocity AMSRE for the 
MLIDDF algorithm is reduced by 23% compared to the UKF and DDF algorithms. The ballistic 
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coefficient AMSRE for the MLIDDF algorithm is reduced by 7% compared to the UKF and DDF 
algorithms. Although the ballistic coefficient AMSRE for the MLIDDF algorithm has no significant 
reduction, we can see that it is reduced significantly compared to the UKF, DDF, lUKF and IDDF 
algorithms, hence, the MLIDDF is preferred over the other filters in the light of position, velocity and 
ballistic coefficient AMSREs. 

From Table 3, we can see the runtime of MLIDDF is less than that of the lUKF algorithm, and is 
more than those of UKF, DDF, and IDDF algorithms, so the accuracy of the MLIDDF algorithm is 
improved at the cost of an increased computational burden. Meanwhile, we can observe the AMSREs 
for the UKF and DDF in position, velocity and ballistic coefficient are almost identical, and the 
AMSREs of the lUKF and IDDF are almost the same. Therefore, on the basis of the simulation results 
presented in Figures 6-8 and Table 2, one can draw the conclusion that the MLIDDF yields a superior 
performance over the other filters. 

4. Conclusions and Future Work 

In this study, we provide the maximum likelihood based iterated divided difference filter which 
inherits the virtues of the divided difference filter and contains the iteration process in the 
measurement update step. The sequence obtained is guaranteed to move up the likelihood surface 
using the iteration termination condition based on the maximum likelihood surface. The maximum 
likelihood based iterated divided difference is implemented easily and is derivative-free. We apply the 
new filter to state estimation for a ballistic target re-entry scenario and compare its performance 
against the unscented Kalman filter, divided difference filter, iterated unscented Kalman filter and 
iterated divided difference filter with the traditional termination criteria. Simulation results 
demonstrate that the maximum likelihood-based iterated divided difference is much more effective 
than the other filters. The maximum likelihood-based iterated divided difference greatly improves the 
performance of state estimation and has a shorter convergence time. 

Future work may focus on the applications of the maximum likelihood iteration divided difference 
filter to remove the outliers which is a serious deviation from the sample and caused by blink and 
subjective eye movement in video nystagmus signal samples of pilot candidates. 
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